home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / pascal / p_mat.exe / PMATOP.PAS < prev    next >
Pascal/Delphi Source File  |  1993-01-17  |  48KB  |  1,555 lines
  1. {
  2. P-Mat - A Turbo Pascal program for Recursive Matrix Algebra
  3.  
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/30/93
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9.  
  10.  
  11. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  12. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  13. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  14. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  15. FROM USE OF THIS PROGRAM.
  16. }
  17.  
  18. Unit pmatop;
  19.  
  20. Interface
  21. Uses pmat;
  22.  
  23. Function MSort( a: vmatrixptr; col, order: integer ): vmatrixptr;
  24. Function Mexp( ROp: vmatrixptr ): vmatrixptr;
  25. Function Mabs( ROp: vmatrixptr ): vmatrixptr;
  26. Function Mlog( ROp: vmatrixptr ): vmatrixptr;
  27. Function Msqrt( ROp: vmatrixptr ): vmatrixptr;
  28.  
  29. Function Trace( ROp: vmatrixptr ): double;
  30. Function Helm( n: integer ): vmatrixptr;
  31. Function Index( lower, upper, step: integer ): vmatrixptr;
  32. Function Kron( a, b: vmatrixptr ): vmatrixptr;
  33. Function Det( Datain: vmatrixptr ): double;
  34.  
  35. Function Cholesky( ROp: vmatrixptr ): vmatrixptr;
  36. Procedure QR( Var ROp, Q, R: vmatrixptr; makeq: boolean );
  37. Function Svd( Var A, U, S, V: vmatrixptr; 
  38.              makeu, makev: boolean ) : integer;
  39. Function Ginv( ROp: vmatrixptr ): vmatrixptr;
  40. Function Fft( ROp: vmatrixptr; INZEE: boolean ): vmatrixptr;
  41.  
  42. Function Vec( ROp: vmatrixptr ): vmatrixptr;
  43. Function Diag( ROp: vmatrixptr ): vmatrixptr;
  44. Function Shape( ROp: vmatrixptr; nrows: integer ): vmatrixptr;
  45.  
  46. Type  margins = (ALL,ROWS,COLUMNS);
  47. Function Sum( ROp: vmatrixptr; marg: margins ) : vmatrixptr;
  48. Function Sumsq( ROp: vmatrixptr; marg: margins ): vmatrixptr;
  49. Function Cusum( ROp: vmatrixptr ): vmatrixptr;
  50. Function Mmax( ROp: vmatrixptr; marg: margins ): vmatrixptr;
  51. Function Mmin( ROp: vmatrixptr; marg: margins ): vmatrixptr;
  52.  
  53. Procedure Crow( Var ROp: vmatrixptr; row: integer; c: double );
  54. Procedure Srow( Var ROp: vmatrixptr; row1, row2: integer );
  55. Procedure Lrow( Var ROp: vmatrixptr; row1, row2: integer; c: double );
  56. Procedure Ccol( Var ROp: vmatrixptr; col: integer; c: double );
  57. Procedure Scol( Var ROp: vmatrixptr; col1, col2: integer );
  58. Procedure Lcol( Var ROp: vmatrixptr; col1, col2: integer; c: double );
  59.  
  60.  
  61. (************************************************************)
  62. Implementation
  63.  
  64. Procedure heapsort( Var a: vmatrixptr );
  65. Var
  66.   n, s0, s1, s2, s3, s4, s5, s4m1, s5m1 : integer;
  67.   t1, ts : double;
  68. Begin
  69.     { Shell-Metzger sort, PC-World, May 1986 }
  70.     n := a^.r;
  71.     s0 := n;
  72.     s1 := n;
  73.     s1 := s1 Div 2;
  74.     While s1 > 0 Do Begin
  75.         s2 := s0 - s1;
  76.         For s3 := 1 To s2 Do Begin 
  77.             s4 := s3;
  78.             While s4 > 0 Do Begin 
  79.                 s5 := s4 + s1;
  80.                 s4m1 := s4;
  81.                 s5m1 := s5;
  82.                 If a^.m( s4m1, 1 ) > a^.m( s5m1, 1 ) Then Begin 
  83.                     t1 := a^.m( s4m1, 1 );
  84.                     a^.mm( s4m1, 1 )^ := a^.m( s5m1, 1 );
  85.                     a^.mm( s5m1, 1 )^ := t1;
  86.                     ts := a^.m( s4m1, 2 );
  87.                     a^.mm( s4m1, 2 )^ := a^.m( s5m1, 2 );
  88.                     a^.mm( s5m1, 2 )^ := ts;
  89.                     s4 := s4 - s1;
  90.                 End
  91.                 Else Begin 
  92.                     s4 := 0;
  93.                 End;
  94.             End;
  95.         End;
  96.         s1 := s1 Div 2;
  97.     End;
  98. End; { end heapsort }
  99.  
  100.  
  101. Function MSort{(a :vmatrixptr; col, order: integer): vmatrixptr;};
  102. Var
  103.   sorted, temp : vmatrixptr;
  104.   i,j,t : integer;
  105. Begin
  106.     a^.Garbage;
  107.     new( sorted, makematrix( a^.r, a^.c ) );
  108.     If order <> 0 Then Begin 
  109.         { sort column col, then reorder other rows }
  110.         new( temp, makematrix( a^.r, 2 ) );
  111.         For i := 1 To a^.r Do Begin 
  112.             temp^.mm( i, 1 )^ := a^.m( i, col );
  113.             temp^.mm( i, 2 )^ := i ;
  114.         End;
  115.         heapsort( temp );
  116.         For i := 1 To a^.r Do Begin 
  117.             For j := 1 To a^.c Do Begin 
  118.                 t := Trunc( temp^.m( i, 2 ) );
  119.                 sorted^.mm( i, j )^ := a^.m( t, j );
  120.             End;
  121.         End;
  122.     End
  123.     Else Begin 
  124.         { sort row col, then reorder other columns }
  125.         new( temp, makematrix( a^.c, 2 ) );
  126.         For i := 1 To a^.c Do Begin 
  127.             temp^.mm( i, 1 )^ := a^.m( col, i );
  128.             temp^.mm( i, 2 )^ := i ;
  129.         End;
  130.         heapsort( temp );
  131.         For i := 1 To a^.c Do Begin 
  132.             For j := 1 To a^.r Do Begin  
  133.                 t := Trunc( temp^.m( i, 2 ) );
  134.                 sorted^.mm( j, i )^ := a^.m( j, t );
  135.             End;
  136.         End;
  137.     End;
  138.     Dispatch^.Push( sorted );
  139.     dispose( temp, killVmatrix );
  140.     MSort := Dispatch^.ReturnMat;
  141. End;
  142. Function Mexp{(ROp :vmatrixptr): vmatrixptr;};
  143. Var
  144.   temp : vmatrixptr;
  145.   i,j  : integer;
  146. Begin
  147.     {  take exponent of all elements  }
  148.     ROp^.Garbage;
  149.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  150.     For i := 1 To  ROp^.r Do Begin 
  151.         For j := 1 To  ROp^.c Do
  152.             temp^.mm( i, j )^ := exp( ROp^.m( i, j ) );
  153.     End;
  154.     Dispatch^.Push( temp );
  155.     Mexp := Dispatch^.ReturnMat;
  156. End;
  157.  
  158. Function Mabs{(ROp: vmatrixptr): vmatrixptr;};
  159. Var
  160.   temp : vmatrixptr;
  161.   i,j  : integer;
  162. Begin
  163.     {  take absolute values of all elements  }
  164.     ROp^.Garbage;
  165.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  166.     For i := 1 To  ROp^.r Do Begin    
  167.         For j := 1 To  ROp^.c Do
  168.             temp^.mm( i, j )^ := abs( ROp^.m( i, j ) );
  169.     End;
  170.     Dispatch^.Push( temp );
  171.     Mabs := Dispatch^.ReturnMat;
  172. End;
  173.  
  174. Function Mlog{(ROp :vmatrixptr): vmatrixptr;};
  175. Var
  176.   temp : vmatrixptr;
  177.   i,j  : integer;
  178.   R    : double;
  179. Begin
  180.     {  take logarithm of all elements  }
  181.     ROp^.Garbage;
  182.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  183.     For i := 1 To  ROp^.r Do Begin    
  184.         For j := 1 To  ROp^.c Do Begin 
  185.             R := ROp^.m( i, j );
  186.             If R <= 0.0 Then
  187.                 Dispatch^.Nrerror( 'Mlog: zero or negative argument to log' )
  188.             Else temp^.mm( i, j )^ := ln( R );
  189.         End;
  190.     End;
  191.     Dispatch^.Push( temp );
  192.     Mlog := Dispatch^.ReturnMat;
  193. End;
  194.  
  195. Function Msqrt{(ROp: vmatrixptr): vmatrixptr;};
  196. Var
  197.   temp : vmatrixptr;
  198.   i,j  : integer;
  199.   R    : double;
  200. Begin
  201.     {  take sqrt of all elements  }
  202.     ROp^.Garbage;
  203.     new( temp, makematrix( ROp^.r, ROp^.c ) );
  204.     For i := 1 To  ROp^.r Do Begin 
  205.         For j := 1 To  ROp^.c Do Begin 
  206.             R := ROp^.m( i, j );
  207.             If R < 0 Then
  208.                 Dispatch^.Nrerror( 'Msqrt: zero or negative argument to sqrt' )
  209.             Else temp^.mm( i, j )^ := sqrt( R );
  210.         End;
  211.     End;
  212.     Dispatch^.Push( temp );
  213.     Msqrt := Dispatch^.ReturnMat;
  214. End;
  215.  
  216. Function Trace{(ROp: vmatrixptr): double;};
  217. Var
  218.   i : integer;
  219.   t : double;
  220. Begin
  221.     ROp^.Garbage;
  222.     t := 0;
  223.     If ROp^.r <> ROp^.c Then
  224.         Dispatch^.Nrerror( ' Trace: matrix not square in trace' );
  225.     For i := 1 To  ROp^.r Do t := t + ROp^.m( i, i );
  226.     trace := t;
  227. End;
  228.  
  229. Function Helm {(n: integer): vmatrixptr;};
  230. Var
  231.   i, j : integer;
  232.   d, den : double;
  233.   temp : vmatrixptr;
  234. Begin
  235.     new( temp, makematrix( n, n ) );
  236.     For i := 1  To n Do Begin    
  237.         For j := 1 To  n Do Begin    
  238.             d := 1.0 / sqrt( 0.0 + n );
  239.             den := d;
  240.             If j > 1 Then den := 1.0 / sqrt( 0.0 + j * (j - 1) );
  241.             If (j > 1) And  (j < i) Then d := 0;
  242.             If (j > 1) And  (j = i) Then d := - den * ( 0.0 + (j - 1));
  243.             If (j > 1) And  (j > i) Then d := den;
  244.             temp^.mm( i, j )^ := d;
  245.         End;
  246.     End;
  247.     Dispatch^.Push( temp );
  248.     helm := Dispatch^.ReturnMat;
  249. End;
  250.  
  251. Function Index{( lower,  upper, step : integer): vmatrixptr;};
  252. Var
  253.   cnter, i : integer;
  254.   temp : vmatrixptr;
  255. Begin
  256.     If step = 0 Then Dispatch^.Nrerror( 'Index: step is zero' );
  257.     cnter := 0;
  258.     i := lower;
  259.     If lower < upper Then Begin 
  260.         If step < 0 Then
  261.             Dispatch^.Nrerror( 'Index: trying to step from lower to upper with negative step size' );
  262.         While i <= upper Do Begin 
  263.             cnter := cnter + 1;
  264.             i := i + step;
  265.         End;
  266.     End
  267.     Else Begin 
  268.